We present a rice-wheat systen spatial decisionb support tool which relies on data from the matched landscape crop assessment surveys in eastern India and a multivariate geoadditive Bayesian model to predict entry points for system optimization in the area of interest.
library(sp)
library(mgcv)
library(bamlss)
# Multivariate geoadditive model
# remotes::install_git("https://git.uibk.ac.at/c4031039/mvnchol")
# library(mvnchol)
library(BayesX)
library(R2BayesX)
library(sf)
library(spdep)
library(rio)
Irrig_Rev_rice_wheat <- import("data/Irrig_Rev_rice_wheat.csv")
library(janitor)
# Irrig_Rev_rice_wheat <- clean_names(Irrig_Rev_rice_wheat)
library(tidyr)
#Irrig_Rev_rice_wheat <- Irrig_Rev_rice_wheat %>% drop_na()
# library(lubridate)
# library(anytime)
# Irrig_Rev_rice_wheat$januaryfirst2017 <- ymd("2017-01-01")
# Irrig_Rev_rice_wheat$sowdate_fmt_wheat_day <- Irrig_Rev_rice_wheat$sowdate_fmt - Irrig_Rev_rice_wheat$januaryfirst2017
# Irrig_Rev_rice_wheat$sowdate_fmt_rice_day <- Irrig_Rev_rice_wheat$sowdate_fmt_rice - Irrig_Rev_rice_wheat$januaryfirst2017
Irrig_Rev_rice_wheat$harvest_day_rice <- Irrig_Rev_rice_wheat$l_crop_duration_days_rice + Irrig_Rev_rice_wheat$sowdate_fmt_rice_day
# Irrig_Rev_rice_wheat$harvest_day_wheat <- Irrig_Rev_rice_wheat$l_crop_duration_days + Irrig_Rev_rice_wheat$sowdate_fmt_wheat_day
shpname <- file.path(getwd(), "shp", "India_aoi_sf_sp")
India_aoi_sp_bnd <- BayesX::shp2bnd(shpname = shpname, regionnames = "District", check.is.in = F)Reading map ... finished
Note: map consists originally of 50 polygons
Note: map consists of 47 regions
f_rice_yield_MRF <- list(
b_grain_yield_ton_per_ha_rice ~ 1 + rice_duration_class_long + s(sowdate_fmt_rice_day) + s(g_q5305_irrig_times_rice) + s(nperha_rice) + s(p2o5perha_rice) + s(District, bs = "mrf", xt = list("penalty" = K)) +
s(District, bs = "re"),
sigma ~ 1 + rice_duration_class_long + s(sowdate_fmt_rice_day) + s(g_q5305_irrig_times_rice) + s(nperha_rice) + s(p2o5perha_rice) + s(District, bs = "mrf", xt = list("penalty" = K)) +
s(District, bs = "re")
)
K <- neighbormatrix(India_aoi_sp_bnd)
head(K) Araria Arwal Aurangabad Banka Begusarai Bhagalpur Arah Buxar
Araria 4 0 0 0 0 0 0 0
Arwal 0 6 -1 0 0 0 -1 0
Aurangabad 0 -1 3 0 0 0 0 0
Banka 0 0 0 3 0 -1 0 0
Begusarai 0 0 0 0 5 0 0 0
Bhagalpur 0 0 0 -1 0 6 0 0
Darbhanga Gaya Gopalganj Jamui Jehanabad Kaimur Katihar Khagaria
Araria 0 0 0 0 0 0 0 0
Arwal 0 -1 0 0 -1 0 0 0
Aurangabad 0 -1 0 0 0 0 0 0
Banka 0 0 0 -1 0 0 0 0
Begusarai 0 0 0 0 0 0 0 -1
Bhagalpur 0 0 0 0 0 0 -1 -1
Kishanganj Lakhisarai Madhepura Madhubani Munger Muzaffarpur Nalanda
Araria -1 0 -1 0 0 0 0
Arwal 0 0 0 0 0 0 0
Aurangabad 0 0 0 0 0 0 0
Banka 0 0 0 0 -1 0 0
Begusarai 0 -1 0 0 -1 0 0
Bhagalpur 0 0 -1 0 -1 0 0
Nawada WestChamparan Patna EastChamparan Purnia Rohtas Saharsa
Araria 0 0 0 0 -1 0 0
Arwal 0 0 -1 0 0 -1 0
Aurangabad 0 0 0 0 0 -1 0
Banka 0 0 0 0 0 0 0
Begusarai 0 0 -1 0 0 0 0
Bhagalpur 0 0 0 0 -1 0 0
Samastipur Saran Sheikhpura Sheohar Sitamarhi Siwan Supaul Vaishali
Araria 0 0 0 0 0 0 -1 0
Arwal 0 0 0 0 0 0 0 0
Aurangabad 0 0 0 0 0 0 0 0
Banka 0 0 0 0 0 0 0 0
Begusarai -1 0 0 0 0 0 0 0
Bhagalpur 0 0 0 0 0 0 0 0
Balia Chandauli Deoria Gazipur Gorakhpur Kushinagar Maharajganj Mau
Araria 0 0 0 0 0 0 0 0
Arwal 0 0 0 0 0 0 0 0
Aurangabad 0 0 0 0 0 0 0 0
Banka 0 0 0 0 0 0 0 0
Begusarai 0 0 0 0 0 0 0 0
Bhagalpur 0 0 0 0 0 0 0 0
Siddharthnagar
Araria 0
Arwal 0
Aurangabad 0
Banka 0
Begusarai 0
Bhagalpur 0
## Also need to transform to factor for
## setting up the MRF smooth.
Irrig_Rev_rice_wheat$District <- as.factor(Irrig_Rev_rice_wheat$a_q103_district)
## Now note that not all regions are observed,
## therefore we need to remove those regions
## from the penalty matrix
rn <- rownames(K)
lv <- levels(Irrig_Rev_rice_wheat$District)
i <- rn %in% lv
K <- K[i, i]
set.seed(321)
library(bamlss)
b_rice_yield_MRF <- bamlss(f_rice_yield_MRF, data = Irrig_Rev_rice_wheat, family = "gaussian")AICc 12433.92 logPost -6216.75 logLik -6063.08 edf 148.79 eps 2.3399 iteration 1
AICc 12239.92 logPost -5910.29 logLik -5967.69 edf 147.28 eps 0.4613 iteration 2
AICc 12211.69 logPost -5884.01 logLik -5956.93 edf 144.14 eps 0.1461 iteration 3
AICc 12205.11 logPost -5882.70 logLik -5953.35 edf 144.42 eps 2.9473 iteration 4
AICc 12203.12 logPost -5882.46 logLik -5952.28 edf 144.48 eps 0.0429 iteration 5
AICc 12202.50 logPost -5882.43 logLik -5951.95 edf 144.50 eps 0.0113 iteration 6
AICc 12202.32 logPost -5882.45 logLik -5951.84 edf 144.52 eps 0.0088 iteration 7
AICc 12202.28 logPost -5882.44 logLik -5951.82 edf 144.52 eps 0.0010 iteration 8
AICc 12202.27 logPost -5882.44 logLik -5951.81 edf 144.53 eps 0.0002 iteration 9
AICc 12202.27 logPost -5882.44 logLik -5951.81 edf 144.53 eps 0.0001 iteration 10
AICc 12202.27 logPost -5882.44 logLik -5951.81 edf 144.53 eps 0.0001 iteration 11
AICc 12202.26 logPost -5882.44 logLik -5951.80 edf 144.53 eps 0.0000 iteration 12
AICc 12202.26 logPost -5882.44 logLik -5951.80 edf 144.53 eps 0.0000 iteration 12
elapsed time: 7.82sec
Starting the sampler...
| | 0% 1.19min
|* | 5% 1.21min 3.83sec
|** | 10% 1.15min 7.66sec
|*** | 15% 1.10min 11.66sec
|**** | 20% 1.05min 15.77sec
|***** | 25% 1.01min 20.13sec
|****** | 30% 57.10sec 24.47sec
|******* | 35% 53.63sec 28.88sec
|******** | 40% 49.71sec 33.14sec
|********* | 45% 45.76sec 37.44sec
|********** | 50% 41.75sec 41.75sec
|*********** | 55% 37.70sec 46.08sec
|************ | 60% 33.66sec 50.49sec
|************* | 65% 29.48sec 54.75sec
|************** | 70% 25.32sec 59.07sec
|*************** | 75% 21.13sec 1.06min
|**************** | 80% 16.94sec 1.13min
|***************** | 85% 12.73sec 1.20min
|****************** | 90% 8.48sec 1.27min
|******************* | 95% 4.23sec 1.34min
|********************| 100% 0.00sec 1.41min
## First, note that we have the structured id = 'mrf1' and unstructured
## spatial effect id = 're2', also indicated in the model summary
summary(b_rice_yield_MRF)
Call:
bamlss(formula = f_rice_yield_MRF, family = "gaussian", data = Irrig_Rev_rice_wheat)
---
Family: gaussian
Link function: mu = identity, sigma = log
*---
Formula mu:
---
b_grain_yield_ton_per_ha_rice ~ 1 + rice_duration_class_long +
s(sowdate_fmt_rice_day) + s(g_q5305_irrig_times_rice) + s(nperha_rice) +
s(p2o5perha_rice) + s(District, bs = "mrf", xt = list(penalty = K)) +
s(District, bs = "re")
-
Parametric coefficients:
Mean 2.5% 50% 97.5% parameters
(Intercept) 1.90903 0.53256 1.64495 3.75492 0.013
rice_duration_class_long 0.07427 0.01184 0.07438 0.13788 0.075
-
Acceptance probability:
Mean 2.5% 50% 97.5%
alpha 1 1 1 1
-
Smooth terms:
Mean 2.5% 50% 97.5%
s(sowdate_fmt_rice_day).tau21 5.510e-01 1.460e-02 2.985e-01 2.446e+00
s(sowdate_fmt_rice_day).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(sowdate_fmt_rice_day).edf 3.784e+00 1.625e+00 3.700e+00 5.909e+00
s(g_q5305_irrig_times_rice).tau21 2.324e+00 2.653e-01 1.653e+00 8.760e+00
s(g_q5305_irrig_times_rice).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(g_q5305_irrig_times_rice).edf 6.482e+00 4.641e+00 6.530e+00 7.994e+00
s(nperha_rice).tau21 4.816e-01 1.694e-04 7.272e-02 3.602e+00
s(nperha_rice).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(nperha_rice).edf 2.866e+00 1.010e+00 2.509e+00 6.675e+00
s(p2o5perha_rice).tau21 9.498e-02 9.163e-05 7.279e-03 8.451e-01
s(p2o5perha_rice).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(p2o5perha_rice).edf 1.889e+00 1.007e+00 1.411e+00 4.674e+00
s(District,id='mrf1').tau21 4.290e-02 7.425e-05 1.089e-02 2.036e-01
s(District,id='mrf1').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(District,id='mrf1').edf 2.274e+01 1.975e+00 2.845e+01 3.410e+01
s(District,id='re2').tau21 5.664e+00 2.034e-01 5.310e+00 1.506e+01
s(District,id='re2').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(District,id='re2').edf 3.560e+01 3.391e+01 3.581e+01 3.594e+01
parameters
s(sowdate_fmt_rice_day).tau21 0.182
s(sowdate_fmt_rice_day).alpha NA
s(sowdate_fmt_rice_day).edf 3.291
s(g_q5305_irrig_times_rice).tau21 2.895
s(g_q5305_irrig_times_rice).alpha NA
s(g_q5305_irrig_times_rice).edf 7.106
s(nperha_rice).tau21 1.398
s(nperha_rice).alpha NA
s(nperha_rice).edf 5.489
s(p2o5perha_rice).tau21 0.018
s(p2o5perha_rice).alpha NA
s(p2o5perha_rice).edf 1.787
s(District,id='mrf1').tau21 0.029
s(District,id='mrf1').alpha NA
s(District,id='mrf1').edf 31.761
s(District,id='re2').tau21 12.212
s(District,id='re2').alpha NA
s(District,id='re2').edf 35.917
---
Formula sigma:
---
sigma ~ 1 + rice_duration_class_long + s(sowdate_fmt_rice_day) +
s(g_q5305_irrig_times_rice) + s(nperha_rice) + s(p2o5perha_rice) +
s(District, bs = "mrf", xt = list(penalty = K)) + s(District,
bs = "re")
-
Parametric coefficients:
Mean 2.5% 50% 97.5% parameters
(Intercept) -0.119970 -0.179496 -0.119649 -0.063475 -0.114
rice_duration_class_long 0.061685 0.009405 0.061474 0.115939 0.060
-
Acceptance probability:
Mean 2.5% 50% 97.5%
alpha 0.9860 0.8949 1.0000 1
-
Smooth terms:
Mean 2.5% 50% 97.5%
s(sowdate_fmt_rice_day).tau21 2.255e-02 9.588e-05 4.358e-03 1.617e-01
s(sowdate_fmt_rice_day).alpha 9.779e-01 8.326e-01 1.000e+00 1.000e+00
s(sowdate_fmt_rice_day).edf 1.623e+00 1.009e+00 1.340e+00 3.522e+00
s(g_q5305_irrig_times_rice).tau21 2.200e-01 2.578e-04 1.086e-01 1.030e+00
s(g_q5305_irrig_times_rice).alpha 9.450e-01 6.660e-01 9.945e-01 1.000e+00
s(g_q5305_irrig_times_rice).edf 4.016e+00 1.123e+00 4.209e+00 6.532e+00
s(nperha_rice).tau21 1.741e-01 4.189e-04 3.295e-02 1.329e+00
s(nperha_rice).alpha 9.529e-01 6.319e-01 9.979e-01 1.000e+00
s(nperha_rice).edf 2.684e+00 1.043e+00 2.352e+00 6.037e+00
s(p2o5perha_rice).tau21 3.062e-01 1.414e-04 1.220e-02 2.836e+00
s(p2o5perha_rice).alpha 9.579e-01 6.586e-01 9.982e-01 1.000e+00
s(p2o5perha_rice).edf 2.461e+00 1.017e+00 1.808e+00 6.536e+00
s(District,id='mrf1').tau21 2.374e-03 7.354e-05 1.646e-03 7.531e-03
s(District,id='mrf1').alpha 8.160e-01 2.828e-01 9.140e-01 1.000e+00
s(District,id='mrf1').edf 1.912e+01 2.882e+00 2.035e+01 2.921e+01
s(District,id='re2').tau21 2.233e-02 5.515e-03 2.132e-02 4.458e-02
s(District,id='re2').alpha 7.382e-01 1.551e-01 8.258e-01 1.000e+00
s(District,id='re2').edf 2.842e+01 1.981e+01 2.919e+01 3.199e+01
parameters
s(sowdate_fmt_rice_day).tau21 0.014
s(sowdate_fmt_rice_day).alpha NA
s(sowdate_fmt_rice_day).edf 1.802
s(g_q5305_irrig_times_rice).tau21 0.191
s(g_q5305_irrig_times_rice).alpha NA
s(g_q5305_irrig_times_rice).edf 4.768
s(nperha_rice).tau21 0.091
s(nperha_rice).alpha NA
s(nperha_rice).edf 3.139
s(p2o5perha_rice).tau21 2.253
s(p2o5perha_rice).alpha NA
s(p2o5perha_rice).edf 6.288
s(District,id='mrf1').tau21 0.007
s(District,id='mrf1').alpha NA
s(District,id='mrf1').edf 28.648
s(District,id='re2').tau21 0.002
s(District,id='re2').alpha NA
s(District,id='re2').edf 10.540
---
Sampler summary:
-
DIC = 12123.05 logLik = -6011.865 pd = 99.3237
runtime = 85.53
---
Optimizer summary:
-
AICc = 12202.27 edf = 144.5353 logLik = -5951.809
logPost = -5882.444 nobs = 4537 runtime = 7.82
# Plot the nonlinear effect
plot(b_rice_yield_MRF, model = "mu", term = "s(sowdate_fmt_rice_day)")plot(b_rice_yield_MRF, model = "mu", term = "s(g_q5305_irrig_times_rice)")plot(b_rice_yield_MRF, model = "mu", term = "s(nperha_rice)")plot(b_rice_yield_MRF, model = "mu", term = "s(p2o5perha_rice)")f_wheat_yield_MRF <- list(
l_ton_per_hectare ~ 1 + variety_type_NMWV + s(sowdate_fmt_wheat_day) + g_q5305_irrig_times + nperha + p2o5perha + s(District, bs = "mrf", xt = list("penalty" = K)) +
s(District, bs = "re"),
sigma ~ 1 + variety_type_NMWV + s(sowdate_fmt_wheat_day) + g_q5305_irrig_times + nperha + p2o5perha + s(District, bs = "mrf", xt = list("penalty" = K)) +
s(District, bs = "re")
)
set.seed(321)
b_wheat_yield_MRF <- bamlss(f_wheat_yield_MRF, data = Irrig_Rev_rice_wheat, family = "gaussian")AICc 8958.541 logPost -4419.70 logLik -4335.18 edf 139.61 eps 0.3246 iteration 1
AICc 8445.144 logPost -4096.84 logLik -4071.93 edf 145.76 eps 0.1431 iteration 2
AICc 8344.698 logPost -4042.08 logLik -4030.62 edf 137.40 eps 0.0345 iteration 3
AICc 8319.455 logPost -4032.32 logLik -4022.32 edf 133.33 eps 0.0132 iteration 4
AICc 8313.001 logPost -4030.95 logLik -4020.60 edf 131.91 eps 0.0050 iteration 5
AICc 8312.046 logPost -4029.90 logLik -4020.23 edf 131.81 eps 0.0019 iteration 6
AICc 8311.784 logPost -4029.54 logLik -4020.15 edf 131.76 eps 0.0010 iteration 7
AICc 8311.335 logPost -4029.64 logLik -4020.12 edf 131.57 eps 0.0005 iteration 8
AICc 8311.297 logPost -4029.53 logLik -4020.11 edf 131.57 eps 0.0003 iteration 9
AICc 8311.268 logPost -4029.43 logLik -4020.10 edf 131.56 eps 0.0002 iteration 10
AICc 8311.245 logPost -4029.35 logLik -4020.09 edf 131.56 eps 0.0002 iteration 11
AICc 8311.226 logPost -4029.28 logLik -4020.09 edf 131.56 eps 0.0002 iteration 12
AICc 8311.211 logPost -4029.22 logLik -4020.08 edf 131.56 eps 0.0001 iteration 13
AICc 8311.197 logPost -4029.16 logLik -4020.08 edf 131.55 eps 0.0001 iteration 14
AICc 8311.185 logPost -4029.11 logLik -4020.07 edf 131.55 eps 0.0001 iteration 15
AICc 8311.175 logPost -4029.06 logLik -4020.07 edf 131.55 eps 0.0001 iteration 16
AICc 8311.175 logPost -4029.06 logLik -4020.07 edf 131.55 eps 0.0001 iteration 16
elapsed time: 5.11sec
Starting the sampler...
| | 0% 52.36sec
|* | 5% 53.77sec 2.83sec
|** | 10% 50.58sec 5.62sec
|*** | 15% 48.51sec 8.56sec
|**** | 20% 46.48sec 11.62sec
|***** | 25% 44.49sec 14.83sec
|****** | 30% 42.54sec 18.23sec
|******* | 35% 40.10sec 21.59sec
|******** | 40% 38.51sec 25.67sec
|********* | 45% 37.64sec 30.80sec
|********** | 50% 36.25sec 36.25sec
|*********** | 55% 33.61sec 41.08sec
|************ | 60% 30.01sec 45.02sec
|************* | 65% 26.28sec 48.80sec
|************** | 70% 22.53sec 52.56sec
|*************** | 75% 18.74sec 56.23sec
|**************** | 80% 15.02sec 1.00min
|***************** | 85% 11.28sec 1.07min
|****************** | 90% 7.54sec 1.13min
|******************* | 95% 3.78sec 1.20min
|********************| 100% 0.00sec 1.27min
## First, note that we have the structured id = 'mrf1' and unstructured
## spatial effect id = 're2', also indicated in the model summary
summary(b_wheat_yield_MRF)
Call:
bamlss(formula = f_wheat_yield_MRF, family = "gaussian", data = Irrig_Rev_rice_wheat)
---
Family: gaussian
Link function: mu = identity, sigma = log
*---
Formula mu:
---
l_ton_per_hectare ~ 1 + variety_type_NMWV + s(sowdate_fmt_wheat_day) +
g_q5305_irrig_times + nperha + p2o5perha + s(District, bs = "mrf",
xt = list(penalty = K)) + s(District, bs = "re")
-
Parametric coefficients:
Mean 2.5% 50% 97.5% parameters
(Intercept) -0.5740916 -1.1340813 -0.6630343 0.1609182 -1.359
variety_type_NMWV 0.3056829 0.2611116 0.3052558 0.3529448 0.303
g_q5305_irrig_times 0.4153445 0.3895350 0.4153230 0.4439696 0.411
nperha 0.0012336 0.0006668 0.0012288 0.0017957 0.001
p2o5perha 0.0031334 0.0020527 0.0031371 0.0042060 0.003
-
Acceptance probability:
Mean 2.5% 50% 97.5%
alpha 1 1 1 1
-
Smooth terms:
Mean 2.5% 50% 97.5%
s(sowdate_fmt_wheat_day).tau21 4.037e-02 7.694e-05 3.583e-03 3.627e-01
s(sowdate_fmt_wheat_day).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(sowdate_fmt_wheat_day).edf 1.940e+00 1.015e+00 1.491e+00 4.943e+00
s(District,id='mrf1').tau21 3.294e-03 6.895e-05 9.061e-04 1.507e-02
s(District,id='mrf1').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(District,id='mrf1').edf 1.865e+01 3.933e+00 1.853e+01 3.206e+01
s(District,id='re2').tau21 4.793e+00 1.759e+00 4.533e+00 9.645e+00
s(District,id='re2').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(District,id='re2').edf 3.589e+01 3.577e+01 3.590e+01 3.596e+01
parameters
s(sowdate_fmt_wheat_day).tau21 15.401
s(sowdate_fmt_wheat_day).alpha NA
s(sowdate_fmt_wheat_day).edf 8.309
s(District,id='mrf1').tau21 0.003
s(District,id='mrf1').alpha NA
s(District,id='mrf1').edf 26.415
s(District,id='re2').tau21 6.085
s(District,id='re2').alpha NA
s(District,id='re2').edf 35.897
---
Formula sigma:
---
sigma ~ 1 + variety_type_NMWV + s(sowdate_fmt_wheat_day) + g_q5305_irrig_times +
nperha + p2o5perha + s(District, bs = "mrf", xt = list(penalty = K)) +
s(District, bs = "re")
-
Parametric coefficients:
Mean 2.5% 50% 97.5% parameters
(Intercept) -0.7819694 -0.8953612 -0.7812289 -0.6696800 -0.316
variety_type_NMWV 0.2056530 0.1502403 0.2059671 0.2547450 0.198
g_q5305_irrig_times 0.0927109 0.0623800 0.0929975 0.1245421 0.098
nperha -0.0011626 -0.0018926 -0.0011700 -0.0004521 -0.001
p2o5perha 0.0017864 0.0003726 0.0017712 0.0031635 0.002
-
Acceptance probability:
Mean 2.5% 50% 97.5%
alpha 0.8957 0.4855 0.9745 1
-
Smooth terms:
Mean 2.5% 50% 97.5%
s(sowdate_fmt_wheat_day).tau21 8.472e-02 5.798e-04 4.162e-02 4.208e-01
s(sowdate_fmt_wheat_day).alpha 9.647e-01 7.699e-01 9.984e-01 1.000e+00
s(sowdate_fmt_wheat_day).edf 2.574e+00 1.067e+00 2.524e+00 4.535e+00
s(District,id='mrf1').tau21 6.406e-04 5.097e-05 3.717e-04 2.418e-03
s(District,id='mrf1').alpha 8.975e-01 4.521e-01 9.764e-01 1.000e+00
s(District,id='mrf1').edf 1.089e+01 2.084e+00 9.704e+00 2.304e+01
s(District,id='re2').tau21 1.986e-02 9.452e-03 1.872e-02 3.562e-02
s(District,id='re2').alpha 7.360e-01 1.839e-01 8.280e-01 1.000e+00
s(District,id='re2').edf 2.834e+01 2.429e+01 2.855e+01 3.129e+01
parameters
s(sowdate_fmt_wheat_day).tau21 0.035
s(sowdate_fmt_wheat_day).alpha NA
s(sowdate_fmt_wheat_day).edf 2.415
s(District,id='mrf1').tau21 0.001
s(District,id='mrf1').alpha NA
s(District,id='mrf1').edf 13.967
s(District,id='re2').tau21 0.203
s(District,id='re2').alpha NA
s(District,id='re2').edf 34.551
---
Sampler summary:
-
DIC = 8219.792 logLik = -4070.916 pd = 77.9595
runtime = 76.61
---
Optimizer summary:
-
AICc = 8311.175 edf = 131.5535 logLik = -4020.075
logPost = -4029.064 nobs = 4537 runtime = 5.11
# Plot the nonlinear effect
plot(b_wheat_yield_MRF, model = "mu", term = "s(sowdate_fmt_wheat_day)")## Now, to predict the spatial effects we set up new data.
nd <- data.frame("District" = unique(Irrig_Rev_rice_wheat$District))
## Predict for the structured spatial effects.
p_str_wheat_yield_MRF <- predict(b_wheat_yield_MRF, newdata = nd, term = "s(District,id='mrf1')", intercept = FALSE)
## And the unstructured spatial effect.
p_unstr_wheat_yield_MRF <- predict(b_wheat_yield_MRF, newdata = nd, term = "s(District,id='re2')", intercept = FALSE)
## MRF smooth effect.
plotmap(India_aoi_sp_bnd,
x = p_str_wheat_yield_MRF$mu, id = nd$District,
main = expression(mu), title = "Structured spatial effect"
)plotmap(India_aoi_sp_bnd,
x = p_str_wheat_yield_MRF$sigma, id = nd$District,
main = expression(sigma), title = "Structured spatial effect"
)## Random effects.
plotmap(India_aoi_sp_bnd,
x = p_unstr_wheat_yield_MRF$mu, id = nd$District, title = "Unstructured spatial effect"
)plotmap(India_aoi_sp_bnd,
x = p_unstr_wheat_yield_MRF$sigma, id = nd$District, title = "Unstructured spatial effect"
)library(sp)
library(mgcv)
library(bamlss)
# Multivariate geoadditive model
# remotes::install_git("https://git.uibk.ac.at/c4031039/mvnchol")
# install_github("https://github.com/meteosimon/mvnchol")
library(mvnchol)
# Remove NAs in the monsoon variables
Irrig_Rev_rice_wheat <- subset(Irrig_Rev_rice_wheat, !(is.na(Irrig_Rev_rice_wheat$onset_2017)))
Irrig_Rev_rice_wheat <- subset(Irrig_Rev_rice_wheat, !(is.na(Irrig_Rev_rice_wheat$monsoon_onset_dev)))
Irrig_Rev_rice_wheat <- subset(Irrig_Rev_rice_wheat, !(is.na(Irrig_Rev_rice_wheat$median_onset_82_15)))
Irrig_Rev_rice_wheat <- subset(Irrig_Rev_rice_wheat, !(is.na(Irrig_Rev_rice_wheat$sd_onset_82_15)))
# Rice first stage
f_sow_rice_1st_stage <- list(
sowdate_fmt_rice_day ~ 1 + rice_duration_class_long + s(gw_2017) + s(onset_2017) + s(monsoon_onset_dev) + s(median_onset_82_15) + s(sd_onset_82_15) + s(District, bs = "mrf", xt = list("penalty" = K)) +
s(District, bs = "re")
)
f_sow_rice_1st_stage_MRF <- bamlss(f_sow_rice_1st_stage, data = Irrig_Rev_rice_wheat, family = "gaussian")AICc 38862.03 logPost -165306. logLik -19371.2 edf 59.010 eps 0.1226 iteration 1
AICc 34841.08 logPost -31072.7 logLik -17339.3 edf 79.788 eps 0.0980 iteration 2
AICc 33498.94 logPost -18399.0 logLik -16648.2 edf 98.970 eps 0.0632 iteration 3
AICc 33180.98 logPost -17093.4 logLik -16491.0 edf 97.287 eps 0.0297 iteration 4
AICc 33096.83 logPost -16986.5 logLik -16452.6 edf 93.751 eps 0.0069 iteration 5
AICc 33056.53 logPost -16935.2 logLik -16435.7 edf 90.674 eps 0.0017 iteration 6
AICc 33040.37 logPost -16899.0 logLik -16429.8 edf 88.587 eps 0.0009 iteration 7
AICc 33035.63 logPost -16869.0 logLik -16428.1 edf 87.868 eps 0.0003 iteration 8
AICc 33028.91 logPost -16845.9 logLik -16427.6 edf 85.110 eps 0.0001 iteration 9
AICc 33028.60 logPost -16840.2 logLik -16427.5 edf 85.077 eps 0.0000 iteration 10
AICc 33028.60 logPost -16840.2 logLik -16427.5 edf 85.077 eps 0.0000 iteration 10
elapsed time: 4.89sec
Starting the sampler...
| | 0% 1.39min
|* | 5% 1.36min 4.31sec
|** | 10% 1.18min 7.89sec
|*** | 15% 1.15min 12.15sec
|**** | 20% 1.10min 16.47sec
|***** | 25% 1.03min 20.58sec
|****** | 30% 57.38sec 24.59sec
|******* | 35% 53.28sec 28.69sec
|******** | 40% 48.93sec 32.62sec
|********* | 45% 43.19sec 35.34sec
|********** | 50% 38.26sec 38.26sec
|*********** | 55% 33.59sec 41.06sec
|************ | 60% 30.39sec 45.58sec
|************* | 65% 26.62sec 49.44sec
|************** | 70% 22.64sec 52.83sec
|*************** | 75% 19.05sec 57.14sec
|**************** | 80% 15.28sec 1.02min
|***************** | 85% 11.52sec 1.09min
|****************** | 90% 7.65sec 1.15min
|******************* | 95% 3.81sec 1.21min
|********************| 100% 0.00sec 1.26min
fitted_f_sow_rice_1st_stage_MRF <- f_sow_rice_1st_stage_MRF$fitted
Irrig_Rev_rice_wheat$Res_rice_sow <- Irrig_Rev_rice_wheat$sowdate_fmt_rice_day - fitted_f_sow_rice_1st_stage_MRF$mu
# Wheat first stage
f_sow_wheat_1st_stage <- list(
sowdate_fmt_wheat_day ~ 1 + variety_type_NMWV + s(harvest_day_rice) + s(gw_2018) + s(District, bs = "mrf", xt = list("penalty" = K)) +
s(District, bs = "re")
)
f_sow_wheat_1st_stage_MRF <- bamlss(f_sow_wheat_1st_stage, data = Irrig_Rev_rice_wheat, family = "gaussian")AICc 49188.65 logPost -488850. logLik -24543.4 edf 50.285 eps 0.2566 iteration 1
AICc 44824.42 logPost -64385.9 logLik -22386.0 edf 26.046 eps 0.0753 iteration 2
AICc 40409.05 logPost -24539.6 logLik -20136.0 edf 67.444 eps 0.0824 iteration 3
AICc 36648.46 logPost -18973.1 logLik -18246.6 edf 76.219 eps 0.0789 iteration 4
AICc 33900.77 logPost -17313.2 logLik -16871.9 edf 77.104 eps 0.0796 iteration 5
AICc 32636.01 logPost -16736.4 logLik -16238.1 edf 78.490 eps 0.0660 iteration 6
AICc 32419.72 logPost -16629.9 logLik -16129.2 edf 79.163 eps 0.0324 iteration 7
AICc 32413.31 logPost -16627.4 logLik -16125.7 edf 79.502 eps 0.0056 iteration 8
AICc 32412.92 logPost -16627.3 logLik -16125.4 edf 79.574 eps 0.0001 iteration 9
AICc 32412.79 logPost -16627.2 logLik -16125.3 edf 79.581 eps 0.0000 iteration 10
AICc 32412.79 logPost -16627.2 logLik -16125.3 edf 79.581 eps 0.0000 iteration 10
elapsed time: 2.08sec
Starting the sampler...
| | 0% 38.08sec
|* | 5% 54.72sec 2.88sec
|** | 10% 59.85sec 6.65sec
|*** | 15% 53.95sec 9.52sec
|**** | 20% 49.20sec 12.30sec
|***** | 25% 44.79sec 14.93sec
|****** | 30% 41.42sec 17.75sec
|******* | 35% 39.24sec 21.13sec
|******** | 40% 36.22sec 24.15sec
|********* | 45% 33.39sec 27.32sec
|********** | 50% 31.19sec 31.19sec
|*********** | 55% 28.46sec 34.79sec
|************ | 60% 25.27sec 37.91sec
|************* | 65% 21.87sec 40.61sec
|************** | 70% 18.45sec 43.05sec
|*************** | 75% 15.01sec 45.04sec
|**************** | 80% 11.77sec 47.08sec
|***************** | 85% 8.68sec 49.16sec
|****************** | 90% 5.70sec 51.27sec
|******************* | 95% 2.82sec 53.58sec
|********************| 100% 0.00sec 55.61sec
fitted_f_sow_wheat_1st_stage_MRF <- f_sow_wheat_1st_stage_MRF$fitted
Irrig_Rev_rice_wheat$Res_wheat_sow <- Irrig_Rev_rice_wheat$sowdate_fmt_wheat_day - fitted_f_sow_wheat_1st_stage_MRF$mu
# Multivariate with sowing dates as endogenous
f_rice_wheat_yield_MRF_corr_endo <- list(
sowdate_fmt_rice_day ~ 1 + rice_duration_class_long + s(gw_2017) + s(onset_2017) + s(monsoon_onset_dev) + s(median_onset_82_15) + s(sd_onset_82_15) + s(District, bs = "mrf", xt = list("penalty" = K)) +
s(District, bs = "re"),
b_grain_yield_ton_per_ha_rice ~ 1 + rice_duration_class_long + s(Res_rice_sow) + s(g_q5305_irrig_times_rice) + s(nperha_rice) + s(p2o5perha_rice) + s(District, bs = "mrf", xt = list("penalty" = K)) +
s(District, bs = "re"),
sowdate_fmt_wheat_day ~ 1 + variety_type_NMWV + s(harvest_day_rice) + s(District, bs = "mrf", xt = list("penalty" = K)) +
s(District, bs = "re"),
l_ton_per_hectare ~ 1 + variety_type_NMWV + s(Res_wheat_sow) + g_q5305_irrig_times + nperha + p2o5perha + s(District, bs = "mrf", xt = list("penalty" = K)) + s(District, bs = "re"),
lamdiag1 ~ 1,
lamdiag2 ~ 1,
lamdiag3 ~ 1,
lamdiag4 ~ 1,
lambda12 ~ 1,
lambda13 ~ 1,
lambda14 ~ 1,
lambda23 ~ 1,
lambda24 ~ s(District, bs = "mrf", xt = list("penalty" = K)) + s(District, bs = "re"),
lambda34 ~ 1
)
multivariate_geo_sow_MRF_corr_endo <- bamlss(f_rice_wheat_yield_MRF_corr_endo, type = "modified", family = mvnchol_bamlss(k = 4), data = Irrig_Rev_rice_wheat)AICc 349476.1 logPost -2543606 logLik -174262. edf 430.14 eps 1.0000 iteration 1
AICc 157743.5 logPost -2447306 logLik -78382.3 edf 441.55 eps 0.1825 iteration 2
AICc 104307.5 logPost -2420173 logLik -51667.3 edf 439.10 eps 0.2070 iteration 3
AICc 89316.74 logPost -2412636 logLik -44189.9 edf 424.43 eps 0.1254 iteration 4
AICc 86452.65 logPost -2411216 logLik -42768.9 edf 415.34 eps 0.0628 iteration 5
AICc 86231.98 logPost -2411101 logLik -42661.9 edf 412.53 eps 0.0202 iteration 6
AICc 86215.73 logPost -2411097 logLik -42659.1 edf 408.13 eps 0.0027 iteration 7
AICc 86213.14 logPost -2411094 logLik -42658.2 edf 407.83 eps 0.0007 iteration 8
AICc 86211.64 logPost -2411092 logLik -42657.7 edf 407.62 eps 0.0003 iteration 9
AICc 86210.65 logPost -2411090 logLik -42657.4 edf 407.43 eps 0.0002 iteration 10
AICc 86205.44 logPost -2411208 logLik -42657.3 edf 405.37 eps 0.0001 iteration 11
AICc 86205.08 logPost -2411207 logLik -42657.3 edf 405.26 eps 0.0001 iteration 12
AICc 86202.45 logPost -2412410 logLik -42657.2 edf 404.20 eps 0.0001 iteration 13
AICc 86202.37 logPost -2412333 logLik -42657.2 edf 404.18 eps 0.0001 iteration 14
AICc 86202.37 logPost -2412333 logLik -42657.2 edf 404.18 eps 0.0001 iteration 14
elapsed time: 55.35sec
Starting the sampler...
| | 0% 8.03min
|* | 5% 7.75min 24.47sec
|** | 10% 7.15min 47.69sec
|*** | 15% 6.83min 1.21min
|**** | 20% 6.57min 1.64min
|***** | 25% 6.39min 2.13min
|****** | 30% 6.05min 2.59min
|******* | 35% 5.64min 3.04min
|******** | 40% 5.27min 3.51min
|********* | 45% 4.84min 3.96min
|********** | 50% 4.42min 4.42min
|*********** | 55% 4.01min 4.90min
|************ | 60% 3.66min 5.49min
|************* | 65% 3.21min 5.96min
|************** | 70% 2.79min 6.51min
|*************** | 75% 2.35min 7.04min
|**************** | 80% 1.87min 7.50min
|***************** | 85% 1.40min 7.95min
|****************** | 90% 56.05sec 8.41min
|******************* | 95% 28.00sec 8.87min
|********************| 100% 0.00sec 9.32min
summary(multivariate_geo_sow_MRF_corr_endo)
Call:
bamlss(formula = f_rice_wheat_yield_MRF_corr_endo, family = mvnchol_bamlss(k = 4),
data = Irrig_Rev_rice_wheat, type = "modified")
---
Family: mvnchol
Link function: mu1 = identity, mu2 = identity, mu3 = identity, mu4 = identity, lamdiag1 = log, lamdiag2 = log, lamdiag3 = log, lamdiag4 = log, lambda12 = identity, lambda13 = identity, lambda14 = identity, lambda23 = identity, lambda24 = identity, lambda34 = identity
*---
Formula mu1:
---
sowdate_fmt_rice_day ~ 1 + rice_duration_class_long + s(gw_2017) +
s(onset_2017) + s(monsoon_onset_dev) + s(median_onset_82_15) +
s(sd_onset_82_15) + s(District, bs = "mrf", xt = list(penalty = K)) +
s(District, bs = "re")
-
Parametric coefficients:
Mean 2.5% 50% 97.5% parameters
(Intercept) 68.6313 65.0608 69.0601 70.4619 2.030
rice_duration_class_long -1.3633 -2.0104 -1.3634 -0.7134 -1.341
-
Acceptance probability:
Mean 2.5% 50% 97.5%
alpha 1 1 1 1
-
Smooth terms:
Mean 2.5% 50% 97.5% parameters
s(gw_2017).tau21 2.095e+00 1.260e-04 1.210e-02 2.125e+01 428.875
s(gw_2017).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00 NA
s(gw_2017).edf 1.326e+00 9.987e-01 1.008e+00 3.388e+00 7.459
s(onset_2017).tau21 7.779e+00 1.169e-04 2.620e-01 7.322e+01 3370.678
s(onset_2017).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00 NA
s(onset_2017).edf 1.790e+00 9.986e-01 1.275e+00 4.686e+00 8.245
s(monsoon_onset_dev).tau21 7.665e+00 2.474e-04 1.551e-01 7.089e+01 16.446
s(monsoon_onset_dev).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00 NA
s(monsoon_onset_dev).edf 1.930e+00 9.996e-01 1.214e+00 5.147e+00 3.687
s(median_onset_82_15).tau21 2.959e+01 8.439e-04 1.024e+01 1.830e+02 18.779
s(median_onset_82_15).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00 NA
s(median_onset_82_15).edf 3.182e+00 1.001e+00 3.276e+00 6.208e+00 3.766
s(sd_onset_82_15).tau21 1.316e+01 1.644e-03 1.563e+00 8.969e+01 0.537
s(sd_onset_82_15).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00 NA
s(sd_onset_82_15).edf 2.170e+00 1.002e+00 1.801e+00 5.006e+00 1.379
s(District,id='mrf1').tau21 3.089e+01 1.575e+01 2.914e+01 5.586e+01 1484.368
s(District,id='mrf1').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00 NA
s(District,id='mrf1').edf 3.430e+01 3.394e+01 3.431e+01 3.457e+01 34.979
s(District,id='re2').tau21 1.606e+04 1.019e+04 1.539e+04 2.517e+04 1.087
s(District,id='re2').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00 NA
s(District,id='re2').edf 3.599e+01 3.599e+01 3.599e+01 3.600e+01 35.240
---
Formula mu2:
---
b_grain_yield_ton_per_ha_rice ~ 1 + rice_duration_class_long +
s(Res_rice_sow) + s(g_q5305_irrig_times_rice) + s(nperha_rice) +
s(p2o5perha_rice) + s(District, bs = "mrf", xt = list(penalty = K)) +
s(District, bs = "re")
-
Parametric coefficients:
Mean 2.5% 50% 97.5% parameters
(Intercept) 2.16281 0.43171 1.97757 3.90909 0.000
rice_duration_class_long 0.10345 0.03601 0.10329 0.17645 0.096
-
Acceptance probability:
Mean 2.5% 50% 97.5%
alpha 1 1 1 1
-
Smooth terms:
Mean 2.5% 50% 97.5%
s(Res_rice_sow).tau21 2.939e+00 7.670e-02 1.377e+00 1.492e+01
s(Res_rice_sow).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(Res_rice_sow).edf 5.543e+00 2.630e+00 5.592e+00 8.087e+00
s(g_q5305_irrig_times_rice).tau21 1.715e+00 2.310e-01 1.227e+00 5.434e+00
s(g_q5305_irrig_times_rice).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(g_q5305_irrig_times_rice).edf 6.171e+00 4.441e+00 6.195e+00 7.609e+00
s(nperha_rice).tau21 7.856e-01 5.357e-03 3.042e-01 3.674e+00
s(nperha_rice).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(nperha_rice).edf 3.921e+00 1.288e+00 3.786e+00 6.668e+00
s(p2o5perha_rice).tau21 2.061e-01 1.334e-04 3.188e-02 1.348e+00
s(p2o5perha_rice).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(p2o5perha_rice).edf 2.348e+00 1.010e+00 2.059e+00 5.205e+00
s(District,id='mrf1').tau21 3.043e-03 8.435e-05 1.067e-03 1.628e-02
s(District,id='mrf1').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(District,id='mrf1').edf 1.462e+01 2.134e+00 1.342e+01 3.017e+01
s(District,id='re2').tau21 4.974e+00 1.229e-01 4.138e+00 1.460e+01
s(District,id='re2').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(District,id='re2').edf 3.529e+01 3.323e+01 3.577e+01 3.593e+01
parameters
s(Res_rice_sow).tau21 16.244
s(Res_rice_sow).alpha NA
s(Res_rice_sow).edf 8.170
s(g_q5305_irrig_times_rice).tau21 1.698
s(g_q5305_irrig_times_rice).alpha NA
s(g_q5305_irrig_times_rice).edf 6.538
s(nperha_rice).tau21 0.991
s(nperha_rice).alpha NA
s(nperha_rice).edf 5.103
s(p2o5perha_rice).tau21 0.172
s(p2o5perha_rice).alpha NA
s(p2o5perha_rice).edf 3.239
s(District,id='mrf1').tau21 0.035
s(District,id='mrf1').alpha NA
s(District,id='mrf1').edf 32.259
s(District,id='re2').tau21 8.948
s(District,id='re2').alpha NA
s(District,id='re2').edf 35.894
---
Formula mu3:
---
sowdate_fmt_wheat_day ~ 1 + variety_type_NMWV + s(harvest_day_rice) +
s(District, bs = "mrf", xt = list(penalty = K)) + s(District,
bs = "re")
-
Parametric coefficients:
Mean 2.5% 50% 97.5% parameters
(Intercept) 114.655 109.483 114.481 120.587 3.830
variety_type_NMWV -3.072 -3.742 -3.074 -2.460 -3.073
-
Acceptance probability:
Mean 2.5% 50% 97.5%
alpha 1 1 1 1
-
Smooth terms:
Mean 2.5% 50% 97.5% parameters
s(harvest_day_rice).tau21 228.435 44.388 160.237 704.149 1035.210
s(harvest_day_rice).alpha 1.000 1.000 1.000 1.000 NA
s(harvest_day_rice).edf 5.934 4.524 5.898 7.396 8.391
s(District,id='mrf1').tau21 54.003 9.448 45.568 154.095 4326.923
s(District,id='mrf1').alpha 1.000 1.000 1.000 1.000 NA
s(District,id='mrf1').edf 34.590 33.755 34.685 34.847 34.994
s(District,id='re2').tau21 49904.116 32148.970 48241.930 80332.511 1.087
s(District,id='re2').alpha 1.000 1.000 1.000 1.000 NA
s(District,id='re2').edf 35.999 35.998 35.999 35.999 35.240
---
Formula mu4:
---
l_ton_per_hectare ~ 1 + variety_type_NMWV + s(Res_wheat_sow) +
g_q5305_irrig_times + nperha + p2o5perha + s(District, bs = "mrf",
xt = list(penalty = K)) + s(District, bs = "re")
-
Parametric coefficients:
Mean 2.5% 50% 97.5% parameters
(Intercept) -0.3451667 -0.9637935 -0.4133652 0.3762429 -1.419
variety_type_NMWV 0.3422512 0.2824899 0.3429785 0.3966410 0.337
g_q5305_irrig_times 0.4219196 0.3932712 0.4221273 0.4494278 0.424
nperha 0.0015496 0.0009547 0.0015550 0.0021134 0.002
p2o5perha 0.0025587 0.0014420 0.0025588 0.0036574 0.002
-
Acceptance probability:
Mean 2.5% 50% 97.5%
alpha 1 1 1 1
-
Smooth terms:
Mean 2.5% 50% 97.5% parameters
s(Res_wheat_sow).tau21 1.861e-02 7.837e-05 1.649e-03 1.452e-01 0.238
s(Res_wheat_sow).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00 NA
s(Res_wheat_sow).edf 1.607e+00 1.011e+00 1.202e+00 3.875e+00 4.795
s(District,id='mrf1').tau21 1.002e-02 5.295e-05 5.635e-03 4.258e-02 0.010
s(District,id='mrf1').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00 NA
s(District,id='mrf1').edf 2.528e+01 2.963e+00 2.938e+01 3.356e+01 31.282
s(District,id='re2').tau21 3.750e+00 1.161e+00 3.640e+00 7.588e+00 5.082
s(District,id='re2').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00 NA
s(District,id='re2').edf 3.586e+01 3.568e+01 3.588e+01 3.594e+01 35.915
---
Formula lamdiag1:
---
lamdiag1 ~ 1
-
Parametric coefficients:
Mean 2.5% 50% 97.5% parameters
(Intercept) -2.216 -2.236 -2.216 -2.194 -2.207
-
Acceptance probability:
Mean 2.5% 50% 97.5%
alpha 0.9896 0.9035 1.0000 1
---
Formula lamdiag2:
---
lamdiag2 ~ 1
-
Parametric coefficients:
Mean 2.5% 50% 97.5% parameters
(Intercept) 0.06651 0.04523 0.06601 0.08788 0.074
-
Acceptance probability:
Mean 2.5% 50% 97.5%
alpha 0.9913 0.9196 1.0000 1
---
Formula lamdiag3:
---
lamdiag3 ~ 1
-
Parametric coefficients:
Mean 2.5% 50% 97.5% parameters
(Intercept) -2.145 -2.167 -2.145 -2.126 -2.141
-
Acceptance probability:
Mean 2.5% 50% 97.5%
alpha 0.9896 0.9163 1.0000 1
---
Formula lamdiag4:
---
lamdiag4 ~ 1
-
Parametric coefficients:
Mean 2.5% 50% 97.5% parameters
(Intercept) 0.5178 0.4983 0.5173 0.5394 0.527
-
Acceptance probability:
Mean 2.5% 50% 97.5%
alpha 0.9910 0.9175 1.0000 1
---
Formula lambda12:
---
lambda12 ~ 1
-
Parametric coefficients:
Mean 2.5% 50% 97.5% parameters
(Intercept) -0.011213 -0.045104 -0.007393 0.028365 -0.004
-
Acceptance probability:
Mean 2.5% 50% 97.5%
alpha 1 1 1 1
---
Formula lambda13:
---
lambda13 ~ 1
-
Parametric coefficients:
Mean 2.5% 50% 97.5% parameters
(Intercept) 0.007309 0.003664 0.007292 0.010696 0.008
-
Acceptance probability:
Mean 2.5% 50% 97.5%
alpha 1 1 1 1
---
Formula lambda14:
---
lambda14 ~ 1
-
Parametric coefficients:
Mean 2.5% 50% 97.5% parameters
(Intercept) 0.01113 0.00202 0.01088 0.01917 0.011
-
Acceptance probability:
Mean 2.5% 50% 97.5%
alpha 1 1 1 1
---
Formula lambda23:
---
lambda23 ~ 1
-
Parametric coefficients:
Mean 2.5% 50% 97.5% parameters
(Intercept) 0.023714 -0.007921 0.023630 0.055818 0.027
-
Acceptance probability:
Mean 2.5% 50% 97.5%
alpha 1 1 1 1
---
Formula lambda24:
---
lambda24 ~ s(District, bs = "mrf", xt = list(penalty = K)) +
s(District, bs = "re")
-
Parametric coefficients:
Mean 2.5% 50% 97.5% parameters
(Intercept) -0.2382 -0.2908 -0.2359 -0.1996 -0.036
-
Acceptance probability:
Mean 2.5% 50% 97.5%
alpha 1 1 1 1
-
Smooth terms:
Mean 2.5% 50% 97.5% parameters
s(District,id='mrf1').tau21 2.387e-03 1.004e-04 2.022e-03 6.647e-03 0.005
s(District,id='mrf1').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00 NA
s(District,id='mrf1').edf 1.479e+01 1.863e+00 1.601e+01 2.430e+01 22.153
s(District,id='re2').tau21 9.553e-03 1.586e-04 6.116e-03 3.340e-02 0.024
s(District,id='re2').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00 NA
s(District,id='re2').edf 1.374e+01 6.443e-01 1.415e+01 2.686e+01 24.454
---
Formula lambda34:
---
lambda34 ~ 1
-
Parametric coefficients:
Mean 2.5% 50% 97.5% parameters
(Intercept) -0.07590 -0.11822 -0.08111 -0.01156 -0.009
-
Acceptance probability:
Mean 2.5% 50% 97.5%
alpha 1 1 1 1
---
Sampler summary:
-
DIC = 85815.46 logLik = -42797.02 pd = 221.4134
runtime = 561.25
---
Optimizer summary:
-
AICc = 86202.37 edf = 404.1805 logLik = -42657.27
logPost = -2412333 nobs = 4527 runtime = 55.35
nd <- data.frame("District" = unique(Irrig_Rev_rice_wheat$District))
# Focusing on the cross-equation correlations
## Predict for the structured spatial effects.
p_str_multivariate_geo_sow_MRF_corr_endo_rice_y <- predict(multivariate_geo_sow_MRF_corr_endo, newdata = nd, term = "s(District,id='mrf1')", intercept = FALSE)
p_str_multivariate_geo_sow_MRF_corr_endo_rice_ydt <- as.data.frame(p_str_multivariate_geo_sow_MRF_corr_endo_rice_y)
## And the unstructured spatial effect.
p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y <- predict(multivariate_geo_sow_MRF_corr_endo, newdata = nd, term = "s(District,id='re2')", intercept = FALSE)
# Rice sowing spatial equation
plotmap(India_aoi_sp_bnd, x = p_str_multivariate_geo_sow_MRF_corr_endo_rice_y$mu1, id = nd$District, main = expression(mu(rice_sowing)), title = "Structured spatial effect")plotmap(India_aoi_sp_bnd, x = p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y$mu1, id = nd$District, main = expression(mu(rice_sowing)), title = "Unstructured spatial effect")# Rice yield spatial equation
plotmap(India_aoi_sp_bnd, x = p_str_multivariate_geo_sow_MRF_corr_endo_rice_y$mu2, id = nd$District, main = expression(mu(rice_yield)), title = "Structured spatial effect")plotmap(India_aoi_sp_bnd, x = p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y$mu2, id = nd$District, main = expression(mu(rice_yield)), title = "Unstructured spatial effect")# Wheat sowing equation
plotmap(India_aoi_sp_bnd, x = p_str_multivariate_geo_sow_MRF_corr_endo_rice_y$mu3, id = nd$District, main = expression(mu(wheat_sowing)), title = "Structured spatial effect")plotmap(India_aoi_sp_bnd, x = p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y$mu3, id = nd$District, main = expression(mu(wheat_sowing)), title = "Unstructured spatial effect")# Wheat yield spatial equation
plotmap(India_aoi_sp_bnd, x = p_str_multivariate_geo_sow_MRF_corr_endo_rice_y$mu4, id = nd$District, main = expression(mu(wheat_yield)), title = "Structured spatial effect")plotmap(India_aoi_sp_bnd, x = p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y$mu4, id = nd$District, main = expression(mu(wheat_yield)), title = "Unstructured spatial effect")# Focusing on the cross-equation correlations
## MRF smooth effect.
plotmap(India_aoi_sp_bnd,
x = p_str_multivariate_geo_sow_MRF_corr_endo_rice_y$lambda24, id = nd$District,
main = expression(lambda(rice, wheat)), title = "Structured spatial effect"
)## Random effects.
plotmap(India_aoi_sp_bnd,
x = p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y$lambda24, id = nd$District, main = expression(lambda(rice, wheat)), title = "Unstructured spatial effect"
)# Rice sowing equation : Non-linear relationships
# s(gw_2017) + s(onset_2017) + s(monsoon_onset_dev) + s(median_onset_82_15) + s(sd_onset_82_15)
plot(multivariate_geo_sow_MRF_corr_endo, model = "mu1", term = "s(gw_2017)")plot(multivariate_geo_sow_MRF_corr_endo, model = "mu1", term = "s(onset_2017)")plot(multivariate_geo_sow_MRF_corr_endo, model = "mu1", term = "s(monsoon_onset_dev)")plot(multivariate_geo_sow_MRF_corr_endo, model = "mu1", term = "s(median_onset_82_15)")plot(multivariate_geo_sow_MRF_corr_endo, model = "mu1", term = "s(sd_onset_82_15)")# Rice yield equation
# s(Res_rice_sow) + s(g_q5305_irrig_times_rice) + s(nperha_rice) + s(p2o5perha_rice)
plot(multivariate_geo_sow_MRF_corr_endo, model = "mu2", term = "s(Res_rice_sow)")plot(multivariate_geo_sow_MRF_corr_endo, model = "mu2", term = "s(g_q5305_irrig_times_rice)")plot(multivariate_geo_sow_MRF_corr_endo, model = "mu2", term = "s(nperha_rice)")plot(multivariate_geo_sow_MRF_corr_endo, model = "mu2", term = "s(p2o5perha_rice)")# Wheat sowing equation
# s(harvest_day_rice) + s(gw_2018)
plot(multivariate_geo_sow_MRF_corr_endo, model = "mu3", term = "s(harvest_day_rice)")# (multivariate_geo_sow_MRF_corr_endo, model = "mu3", term = "s(gw_2018)")
# Wheat yield equation
plot(multivariate_geo_sow_MRF_corr_endo, model = "mu4", term = "s(Res_wheat_sow)")# Fitted values
multivariate_geo_sow_MRF_corr_endo_fitted_values <- multivariate_geo_sow_MRF_corr_endo$fitted
multivariate_geo_sow_MRF_corr_endo_fitted_values <- as.data.frame(multivariate_geo_sow_MRF_corr_endo_fitted_values)
# Merge the fitted results to the data and export
Irrig_Rev_rice_wheat_Mult_Res <- cbind(Irrig_Rev_rice_wheat, multivariate_geo_sow_MRF_corr_endo_fitted_values)
summary(lm(sowdate_fmt_rice_day ~ mu1, Irrig_Rev_rice_wheat_Mult_Res))
Call:
lm(formula = sowdate_fmt_rice_day ~ mu1, data = Irrig_Rev_rice_wheat_Mult_Res)
Residuals:
Min 1Q Median 3Q Max
-44.574 -5.397 -0.269 4.986 44.985
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.03383 3.41385 -0.01 0.992
mu1 1.00018 0.01771 56.48 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9.093 on 4525 degrees of freedom
Multiple R-squared: 0.4135, Adjusted R-squared: 0.4134
F-statistic: 3190 on 1 and 4525 DF, p-value: < 2.2e-16
summary(lm(b_grain_yield_ton_per_ha_rice ~ l_ton_per_hectare, Irrig_Rev_rice_wheat_Mult_Res))
Call:
lm(formula = b_grain_yield_ton_per_ha_rice ~ l_ton_per_hectare,
data = Irrig_Rev_rice_wheat_Mult_Res)
Residuals:
Min 1Q Median 3Q Max
-3.2891 -0.7976 -0.0052 0.7454 3.3894
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.34109 0.05974 55.93 <2e-16 ***
l_ton_per_hectare 0.22766 0.01958 11.63 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.08 on 4525 degrees of freedom
Multiple R-squared: 0.02901, Adjusted R-squared: 0.0288
F-statistic: 135.2 on 1 and 4525 DF, p-value: < 2.2e-16
plot((lm(sowdate_fmt_rice_day ~ mu1, Irrig_Rev_rice_wheat_Mult_Res)))Irrig_Rev_rice_wheat_Mult_Res_sp <- Irrig_Rev_rice_wheat_Mult_Res
coordinates(Irrig_Rev_rice_wheat_Mult_Res_sp) <- c("o_largest_plot_gps_longitude", "o_largest_plot_gps_latitude")
# Map the correlations
# library(modelsummary)
# Irrig_Rev_rice_wheat_Mult_Res_dist <- datasummary(Heading("District") * District ~ Heading("N obs") * N + Heading("%") * Percent() + lambda24 * (Mean + SD), data = Irrig_Rev_rice_wheat_Mult_Res, output = "data.frame")
# library(dplyr)
# Irrig_Rev_rice_wheat_Mult_Res_dist <- rename(Irrig_Rev_rice_wheat_Mult_Res_dist, "Mean_Rice_Wheat_Rho" = "Mean")
# Irrig_Rev_rice_wheat_Mult_Res_dist <- rename(Irrig_Rev_rice_wheat_Mult_Res_dist, "SD_Rice_Wheat_Rho" = "SD")
# Irrig_Rev_rice_wheat_Mult_Res_dist <- subset(Irrig_Rev_rice_wheat_Mult_Res_dist, Irrig_Rev_rice_wheat_Mult_Res_dist$District != "Purnia")
# rice_wheat_yield_rho_dist_sf <- merge(India_aoi_sf, Irrig_Rev_rice_wheat_Mult_Res_dist, by = "District")
# rice_wheat_yield_rho_dist_sf$Mean_Rice_Wheat_Rho <- as.numeric(rice_wheat_yield_rho_dist_sf$Mean_Rice_Wheat_Rho)
# library(mapview)
# mapview(rice_wheat_yield_rho_dist_sf, zcol = "Mean_Rice_Wheat_Rho", layer.name = "Rice wheat equation correlation")
# library(sf)
# rice_wheat_yield_rho_dist_sf_sp <- as_Spatial(rice_wheat_yield_rho_dist_sf)
# library(tmap)
# tmap_mode("view")
# rice_wheat_yield_rho_dist_sf_sp_map <- tm_shape(rice_wheat_yield_rho_dist_sf_sp) +
# tm_polygons(col = "Mean_Rice_Wheat_Rho", title = "Rice wheat equation correlation", style = "quantile") +
# tm_layout(legend.outside = TRUE)
# tmap_save(rice_wheat_yield_rho_dist_sf_sp_map, "figures/rice_wheat_yield_rho_dist_sf_sp_map .png")f_rice_wheat_yield_MRF_corr_endo_General <- list(
sowdate_fmt_rice_day ~ 1 + rice_duration_class_long + s(gw_2017) + s(onset_2017) + s(monsoon_onset_dev) + s(median_onset_82_15) + s(sd_onset_82_15) + s(District, bs = "mrf", xt = list("penalty" = K)) +
s(District, bs = "re"),
b_grain_yield_ton_per_ha_rice ~ 1 + rice_duration_class_long + s(Res_rice_sow) + s(g_q5305_irrig_times_rice) + s(nperha_rice) + s(p2o5perha_rice) + s(District, bs = "mrf", xt = list("penalty" = K)) +
s(District, bs = "re"),
sowdate_fmt_wheat_day ~ 1 + variety_type_NMWV + s(harvest_day_rice) + s(District, bs = "mrf", xt = list("penalty" = K)) +
s(District, bs = "re"),
l_ton_per_hectare ~ 1 + variety_type_NMWV + s(Res_wheat_sow) + g_q5305_irrig_times + nperha + p2o5perha + s(District, bs = "mrf", xt = list("penalty" = K)) + s(District, bs = "re"),
lamdiag1 ~ 1,
lamdiag2 ~ 1,
lamdiag3 ~ 1,
lamdiag4 ~ 1,
lambda12 ~ 1,
lambda13 ~ 1,
lambda14 ~ 1,
lambda23 ~ 1,
lambda24 ~ s(gw_2017)+s(District, bs = "mrf", xt = list("penalty" = K)) + s(District, bs = "re"),
lambda34 ~ 1
)
multivariate_geo_sow_MRF_corr_endo_General <- bamlss(f_rice_wheat_yield_MRF_corr_endo_General, type = "modified", family = mvnchol_bamlss(k = 4), data = Irrig_Rev_rice_wheat)AICc 349479.0 logPost -2543604 logLik -174260. edf 432.85 eps 1.0000 iteration 1
AICc 157745.3 logPost -2447296 logLik -78381.2 edf 443.21 eps 0.1854 iteration 2
AICc 104308.9 logPost -2420161 logLik -51666.0 edf 440.78 eps 0.2158 iteration 3
AICc 89318.46 logPost -2412624 logLik -44188.5 edf 426.28 eps 0.1451 iteration 4
AICc 86454.46 logPost -2411205 logLik -42767.5 edf 417.25 eps 0.0645 iteration 5
AICc 86233.78 logPost -2411090 logLik -42660.5 edf 414.43 eps 0.0206 iteration 6
AICc 86217.50 logPost -2411085 logLik -42657.7 edf 410.02 eps 0.0028 iteration 7
AICc 86214.90 logPost -2411082 logLik -42656.8 edf 409.71 eps 0.0010 iteration 8
AICc 86213.39 logPost -2411080 logLik -42656.3 edf 409.50 eps 0.0007 iteration 9
AICc 86212.40 logPost -2411078 logLik -42656.0 edf 409.31 eps 0.0045 iteration 10
AICc 86207.19 logPost -2411196 logLik -42655.9 edf 407.25 eps 0.0012 iteration 11
AICc 86206.82 logPost -2411195 logLik -42655.9 edf 407.14 eps 0.0003 iteration 12
AICc 86204.18 logPost -2412413 logLik -42655.8 edf 406.07 eps 0.0003 iteration 13
AICc 86204.09 logPost -2412336 logLik -42655.8 edf 406.05 eps 0.0007 iteration 14
AICc 86204.03 logPost -2412276 logLik -42655.8 edf 406.03 eps 0.0002 iteration 15
AICc 86203.96 logPost -2412276 logLik -42655.8 edf 406.01 eps 0.0001 iteration 16
AICc 86203.91 logPost -2412276 logLik -42655.8 edf 405.99 eps 0.0001 iteration 17
AICc 86203.91 logPost -2412276 logLik -42655.8 edf 405.99 eps 0.0001 iteration 17
elapsed time: 1.07min
Starting the sampler...
| | 0% 6.31min
|* | 5% 5.89min 18.61sec
|** | 10% 5.66min 37.71sec
|*** | 15% 5.40min 57.21sec
|**** | 20% 5.20min 1.30min
|***** | 25% 5.01min 1.67min
|****** | 30% 4.78min 2.05min
|******* | 35% 4.49min 2.42min
|******** | 40% 4.18min 2.79min
|********* | 45% 3.86min 3.16min
|********** | 50% 3.53min 3.53min
|*********** | 55% 3.19min 3.90min
|************ | 60% 2.85min 4.27min
|************* | 65% 2.50min 4.64min
|************** | 70% 2.14min 5.00min
|*************** | 75% 1.79min 5.38min
|**************** | 80% 1.44min 5.74min
|***************** | 85% 1.08min 6.13min
|****************** | 90% 43.32sec 6.50min
|******************* | 95% 21.68sec 6.87min
|********************| 100% 0.00sec 7.24min
summary(multivariate_geo_sow_MRF_corr_endo_General)
Call:
bamlss(formula = f_rice_wheat_yield_MRF_corr_endo_General, family = mvnchol_bamlss(k = 4),
data = Irrig_Rev_rice_wheat, type = "modified")
---
Family: mvnchol
Link function: mu1 = identity, mu2 = identity, mu3 = identity, mu4 = identity, lamdiag1 = log, lamdiag2 = log, lamdiag3 = log, lamdiag4 = log, lambda12 = identity, lambda13 = identity, lambda14 = identity, lambda23 = identity, lambda24 = identity, lambda34 = identity
*---
Formula mu1:
---
sowdate_fmt_rice_day ~ 1 + rice_duration_class_long + s(gw_2017) +
s(onset_2017) + s(monsoon_onset_dev) + s(median_onset_82_15) +
s(sd_onset_82_15) + s(District, bs = "mrf", xt = list(penalty = K)) +
s(District, bs = "re")
-
Parametric coefficients:
Mean 2.5% 50% 97.5% parameters
(Intercept) 71.3421 66.6505 71.2306 76.3630 2.031
rice_duration_class_long -1.3725 -2.0659 -1.3734 -0.7002 -1.342
-
Acceptance probability:
Mean 2.5% 50% 97.5%
alpha 1 1 1 1
-
Smooth terms:
Mean 2.5% 50% 97.5% parameters
s(gw_2017).tau21 2.783e+01 3.293e-02 6.432e+00 1.946e+02 428.733
s(gw_2017).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00 NA
s(gw_2017).edf 2.672e+00 1.023e+00 2.494e+00 5.730e+00 7.459
s(onset_2017).tau21 1.876e+00 6.997e-05 1.240e-02 2.299e+01 3363.637
s(onset_2017).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00 NA
s(onset_2017).edf 1.333e+00 9.976e-01 1.016e+00 3.611e+00 8.244
s(monsoon_onset_dev).tau21 7.569e+00 7.212e-04 4.383e-01 5.875e+01 15.486
s(monsoon_onset_dev).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00 NA
s(monsoon_onset_dev).edf 2.046e+00 1.001e+00 1.480e+00 5.042e+00 3.633
s(median_onset_82_15).tau21 1.339e+01 1.459e-04 4.668e-02 1.026e+02 18.784
s(median_onset_82_15).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00 NA
s(median_onset_82_15).edf 2.045e+00 9.990e-01 1.063e+00 5.504e+00 3.766
s(sd_onset_82_15).tau21 3.826e+00 1.450e-04 2.333e-02 4.081e+01 0.554
s(sd_onset_82_15).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00 NA
s(sd_onset_82_15).edf 1.480e+00 9.990e-01 1.023e+00 4.173e+00 1.388
s(District,id='mrf1').tau21 7.390e+01 3.864e+01 7.183e+01 1.273e+02 1481.527
s(District,id='mrf1').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00 NA
s(District,id='mrf1').edf 3.464e+01 3.443e+01 3.465e+01 3.479e+01 34.979
s(District,id='re2').tau21 1.519e+04 9.677e+03 1.455e+04 2.354e+04 1.087
s(District,id='re2').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00 NA
s(District,id='re2').edf 3.599e+01 3.599e+01 3.599e+01 3.600e+01 35.240
---
Formula mu2:
---
b_grain_yield_ton_per_ha_rice ~ 1 + rice_duration_class_long +
s(Res_rice_sow) + s(g_q5305_irrig_times_rice) + s(nperha_rice) +
s(p2o5perha_rice) + s(District, bs = "mrf", xt = list(penalty = K)) +
s(District, bs = "re")
-
Parametric coefficients:
Mean 2.5% 50% 97.5% parameters
(Intercept) 1.33448 0.07265 1.08753 3.51466 0.001
rice_duration_class_long 0.09706 0.03197 0.09722 0.15912 0.095
-
Acceptance probability:
Mean 2.5% 50% 97.5%
alpha 1 1 1 1
-
Smooth terms:
Mean 2.5% 50% 97.5%
s(Res_rice_sow).tau21 2.419e+00 9.364e-02 1.038e+00 1.288e+01
s(Res_rice_sow).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(Res_rice_sow).edf 5.310e+00 2.798e+00 5.241e+00 7.992e+00
s(g_q5305_irrig_times_rice).tau21 1.857e+00 2.474e-01 1.345e+00 6.468e+00
s(g_q5305_irrig_times_rice).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(g_q5305_irrig_times_rice).edf 6.226e+00 4.500e+00 6.274e+00 7.728e+00
s(nperha_rice).tau21 7.115e-01 1.034e-02 3.575e-01 3.765e+00
s(nperha_rice).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(nperha_rice).edf 3.949e+00 1.483e+00 3.963e+00 6.671e+00
s(p2o5perha_rice).tau21 3.420e-01 1.193e-04 3.836e-02 2.565e+00
s(p2o5perha_rice).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(p2o5perha_rice).edf 2.583e+00 1.008e+00 2.155e+00 5.870e+00
s(District,id='mrf1').tau21 7.992e-03 8.542e-05 2.262e-03 4.889e-02
s(District,id='mrf1').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(District,id='mrf1').edf 1.855e+01 2.203e+00 1.895e+01 3.276e+01
s(District,id='re2').tau21 8.530e+00 4.015e-01 8.099e+00 2.101e+01
s(District,id='re2').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(District,id='re2').edf 3.576e+01 3.472e+01 3.588e+01 3.595e+01
parameters
s(Res_rice_sow).tau21 15.978
s(Res_rice_sow).alpha NA
s(Res_rice_sow).edf 8.159
s(g_q5305_irrig_times_rice).tau21 1.718
s(g_q5305_irrig_times_rice).alpha NA
s(g_q5305_irrig_times_rice).edf 6.551
s(nperha_rice).tau21 0.985
s(nperha_rice).alpha NA
s(nperha_rice).edf 5.097
s(p2o5perha_rice).tau21 0.172
s(p2o5perha_rice).alpha NA
s(p2o5perha_rice).edf 3.240
s(District,id='mrf1').tau21 0.035
s(District,id='mrf1').alpha NA
s(District,id='mrf1').edf 32.263
s(District,id='re2').tau21 8.982
s(District,id='re2').alpha NA
s(District,id='re2').edf 35.894
---
Formula mu3:
---
sowdate_fmt_wheat_day ~ 1 + variety_type_NMWV + s(harvest_day_rice) +
s(District, bs = "mrf", xt = list(penalty = K)) + s(District,
bs = "re")
-
Parametric coefficients:
Mean 2.5% 50% 97.5% parameters
(Intercept) 107.729 105.228 107.925 109.056 3.830
variety_type_NMWV -3.071 -3.671 -3.072 -2.456 -3.073
-
Acceptance probability:
Mean 2.5% 50% 97.5%
alpha 1 1 1 1
-
Smooth terms:
Mean 2.5% 50% 97.5% parameters
s(harvest_day_rice).tau21 249.145 45.941 169.557 943.694 1033.860
s(harvest_day_rice).alpha 1.000 1.000 1.000 1.000 NA
s(harvest_day_rice).edf 5.867 4.288 5.819 7.667 8.391
s(District,id='mrf1').tau21 273.044 152.159 257.988 484.946 4322.922
s(District,id='mrf1').alpha 1.000 1.000 1.000 1.000 NA
s(District,id='mrf1').edf 34.913 34.840 34.914 34.977 34.994
s(District,id='re2').tau21 52391.176 32990.950 50597.124 79970.994 1.087
s(District,id='re2').alpha 1.000 1.000 1.000 1.000 NA
s(District,id='re2').edf 35.998 35.997 35.998 35.999 35.240
---
Formula mu4:
---
l_ton_per_hectare ~ 1 + variety_type_NMWV + s(Res_wheat_sow) +
g_q5305_irrig_times + nperha + p2o5perha + s(District, bs = "mrf",
xt = list(penalty = K)) + s(District, bs = "re")
-
Parametric coefficients:
Mean 2.5% 50% 97.5% parameters
(Intercept) -0.7126342 -1.3986803 -0.8445395 0.2712965 -1.419
variety_type_NMWV 0.3424046 0.2887866 0.3427579 0.3931747 0.337
g_q5305_irrig_times 0.4231647 0.3972220 0.4229748 0.4504678 0.424
nperha 0.0015341 0.0009313 0.0015484 0.0021222 0.002
p2o5perha 0.0025301 0.0013402 0.0025688 0.0036608 0.003
-
Acceptance probability:
Mean 2.5% 50% 97.5%
alpha 1 1 1 1
-
Smooth terms:
Mean 2.5% 50% 97.5% parameters
s(Res_wheat_sow).tau21 2.950e-02 6.382e-05 2.565e-03 2.716e-01 0.233
s(Res_wheat_sow).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00 NA
s(Res_wheat_sow).edf 1.740e+00 1.009e+00 1.296e+00 4.531e+00 4.773
s(District,id='mrf1').tau21 2.583e-03 6.088e-05 5.571e-04 2.028e-02 0.010
s(District,id='mrf1').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00 NA
s(District,id='mrf1').edf 1.662e+01 3.360e+00 1.473e+01 3.267e+01 31.294
s(District,id='re2').tau21 5.518e+00 1.334e+00 5.472e+00 1.160e+01 5.107
s(District,id='re2').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00 NA
s(District,id='re2').edf 3.589e+01 3.571e+01 3.592e+01 3.596e+01 35.915
---
Formula lamdiag1:
---
lamdiag1 ~ 1
-
Parametric coefficients:
Mean 2.5% 50% 97.5% parameters
(Intercept) -2.216 -2.237 -2.216 -2.196 -2.207
-
Acceptance probability:
Mean 2.5% 50% 97.5%
alpha 0.9909 0.9239 1.0000 1
---
Formula lamdiag2:
---
lamdiag2 ~ 1
-
Parametric coefficients:
Mean 2.5% 50% 97.5% parameters
(Intercept) 0.06601 0.04535 0.06603 0.08673 0.074
-
Acceptance probability:
Mean 2.5% 50% 97.5%
alpha 0.9907 0.9241 1.0000 1
---
Formula lamdiag3:
---
lamdiag3 ~ 1
-
Parametric coefficients:
Mean 2.5% 50% 97.5% parameters
(Intercept) -2.145 -2.166 -2.145 -2.125 -2.141
-
Acceptance probability:
Mean 2.5% 50% 97.5%
alpha 0.9905 0.9188 1.0000 1
---
Formula lamdiag4:
---
lamdiag4 ~ 1
-
Parametric coefficients:
Mean 2.5% 50% 97.5% parameters
(Intercept) 0.5179 0.4987 0.5176 0.5392 0.527
-
Acceptance probability:
Mean 2.5% 50% 97.5%
alpha 0.9884 0.9012 1.0000 1
---
Formula lambda12:
---
lambda12 ~ 1
-
Parametric coefficients:
Mean 2.5% 50% 97.5% parameters
(Intercept) -0.019167 -0.039310 -0.018658 0.002473 -0.004
-
Acceptance probability:
Mean 2.5% 50% 97.5%
alpha 1 1 1 1
---
Formula lambda13:
---
lambda13 ~ 1
-
Parametric coefficients:
Mean 2.5% 50% 97.5% parameters
(Intercept) 0.006978 0.003265 0.007021 0.010635 0.008
-
Acceptance probability:
Mean 2.5% 50% 97.5%
alpha 1 1 1 1
---
Formula lambda14:
---
lambda14 ~ 1
-
Parametric coefficients:
Mean 2.5% 50% 97.5% parameters
(Intercept) 0.013333 0.008215 0.013348 0.019123 0.011
-
Acceptance probability:
Mean 2.5% 50% 97.5%
alpha 1 1 1 1
---
Formula lambda23:
---
lambda23 ~ 1
-
Parametric coefficients:
Mean 2.5% 50% 97.5% parameters
(Intercept) 0.023332 -0.008003 0.023810 0.053160 0.027
-
Acceptance probability:
Mean 2.5% 50% 97.5%
alpha 1 1 1 1
---
Formula lambda24:
---
lambda24 ~ s(gw_2017) + s(District, bs = "mrf", xt = list(penalty = K)) +
s(District, bs = "re")
-
Parametric coefficients:
Mean 2.5% 50% 97.5% parameters
(Intercept) -0.2375 -0.2912 -0.2359 -0.2011 -0.037
-
Acceptance probability:
Mean 2.5% 50% 97.5%
alpha 1 1 1 1
-
Smooth terms:
Mean 2.5% 50% 97.5% parameters
s(gw_2017).tau21 3.959e-02 8.403e-05 6.952e-03 2.810e-01 0.019
s(gw_2017).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00 NA
s(gw_2017).edf 1.600e+00 1.005e+00 1.309e+00 3.396e+00 1.644
s(District,id='mrf1').tau21 3.633e-03 2.214e-04 3.565e-03 8.218e-03 0.005
s(District,id='mrf1').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00 NA
s(District,id='mrf1').edf 1.832e+01 3.675e+00 1.986e+01 2.537e+01 22.402
s(District,id='re2').tau21 4.618e-03 1.066e-04 1.174e-03 2.527e-02 0.024
s(District,id='re2').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00 NA
s(District,id='re2').edf 7.958e+00 4.446e-01 4.264e+00 2.514e+01 24.433
---
Formula lambda34:
---
lambda34 ~ 1
-
Parametric coefficients:
Mean 2.5% 50% 97.5% parameters
(Intercept) -0.04277 -0.16725 -0.04010 0.03504 -0.009
-
Acceptance probability:
Mean 2.5% 50% 97.5%
alpha 1 1 1 1
---
Sampler summary:
-
DIC = 85820.91 logLik = -42798.93 pd = 223.049
runtime = 435.95
---
Optimizer summary:
-
AICc = 86203.92 edf = 405.9967 logLik = -42655.86
logPost = -2412277 nobs = 4527 runtime = 64.48
library(distreg.vis)
if (interactive()) {
vis()
}